home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / renaisnc / lib.lha / Lib / meshsurf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-04  |  5.8 KB  |  247 lines
  1. /*
  2.   File: meshsurf.c
  3.   Authors: Charles Loop 
  4.   Last Modified:  June 17 1989
  5.   Purpose: Better implementation (faster, uses less memory)
  6.   */
  7. #include <stdio.h>
  8. #include "Math.h"
  9. #include "4D.h"
  10. #include "mesh.h"
  11.  
  12. extern Plane *PlaneEquation();
  13.  
  14.  
  15. /********** DECasteljau stack variables and constants **********/
  16.  
  17. #define MAX_DEGREE 6
  18. #define STACK_SIZE 84  /* choose(MAX_DEGREE + 3, 3) = (MAX_DEGREE+1)*(MAX_DEGREE+2)*(MAX_DEGREE+3)/6 */
  19.  
  20. static PointType4D B[STACK_SIZE];
  21.  
  22. static int degree;
  23. static int top;
  24.  
  25.  
  26. int LOC(h,i,j,k) /* index function */
  27.      int h,i,j,k;
  28. {
  29.   return top - (k + (h-i)*(h-i+1)/2 + h*(h+1)*(h+2)/6);
  30. }
  31.  
  32. /***************************************************************/
  33.  
  34.  
  35. /************ surface tessellation *****************************/
  36.  
  37. static Vertex *Q;
  38.  
  39. static int EDGE_DIVISIONS;
  40. static int NUMBER_SAMPLES;
  41.  
  42. int ord(i,j,k) /* index function */
  43.      int i,j,k;
  44. {
  45.   return NUMBER_SAMPLES - (EDGE_DIVISIONS-i)*(EDGE_DIVISIONS-i+1)/2 - k - 1;
  46. }
  47.  
  48. /***************************************************************/
  49.  
  50.  
  51. DECasteljau(r,s,t)
  52.      double r,s,t;
  53. {
  54.   
  55.   int h, i, j, k;
  56.   int l0, l1, l2, l3;
  57.   
  58.   for (h = degree-1; h >= 0; h--) {
  59.     for (i = 0; i <= h; i++) {
  60.       for (j = 0; j <= (h-i); j++) {
  61.     k = h - i - j;
  62.     
  63.     l0 = LOC(h,i,j,k);
  64.     l1 = LOC(h+1,i+1,j,k);
  65.     l2 = LOC(h+1,i,j+1,k);
  66.     l3 = LOC(h+1,i,j,k+1);
  67.     
  68.     B[l0].P[X] = r*B[l1].P[X] + s*B[l2].P[X] + t*B[l3].P[X];
  69.     B[l0].P[Y] = r*B[l1].P[Y] + s*B[l2].P[Y] + t*B[l3].P[Y];
  70.     B[l0].P[Z] = r*B[l1].P[Z] + s*B[l2].P[Z] + t*B[l3].P[Z];
  71.     B[l0].P[W] = r*B[l1].P[W] + s*B[l2].P[W] + t*B[l3].P[W];
  72.     
  73.       }
  74.     }
  75.   }
  76. }
  77.  
  78. /***************************************************************/
  79.  
  80. Go(node, r1, r2, r3, s1, s2, s3, t1, t2, t3, l)
  81.      Mesh *node;
  82.      unsigned r1, r2, r3, s1, s2, s3, t1, t2, t3;
  83.      int l;
  84. {
  85.   Mesh *submesh;
  86.   int i;
  87.   
  88.   if (l) {
  89.     node->type = SUBMESH;
  90.     
  91.     submesh = (Mesh *)calloc(4, sizeof(Mesh));
  92.     
  93.     Go(submesh  , (s1+t1)/2, (s2+t2)/2, (s3+t3)/2,
  94.        (t1+r1)/2, (t2+r2)/2, (t3+r3)/2,
  95.        (r1+s1)/2, (r2+s2)/2, (r3+s3)/2, l-1);
  96.     
  97.     Go(submesh+1,  r1,        r2,        r3,
  98.        (r1+s1)/2, (r2+s2)/2, (r3+s3)/2,
  99.        (t1+r1)/2, (t2+r2)/2, (t3+r3)/2, l-1);
  100.     
  101.     Go(submesh+2, (r1+s1)/2, (r2+s2)/2, (r3+s3)/2,
  102.        s1,        s2,        s3,
  103.        (s1+t1)/2, (s2+t2)/2, (s3+t3)/2, l-1);
  104.     
  105.     Go(submesh+3, (t1+r1)/2, (t2+r2)/2, (t3+r3)/2,
  106.        (s1+t1)/2, (s2+t2)/2, (s3+t3)/2,
  107.        t1,        t2,        t3,       l-1);
  108.     
  109.     CombindBoxes(node->box, submesh, submesh+1, submesh+2, submesh+3);
  110.     
  111.     for (i = 0; i < 3; i++)
  112.       submesh[i].next = submesh+i+1;
  113.     
  114.     node->sub.m = submesh;
  115.     
  116.   }
  117.   else {
  118.     node->type = TRIANGLE;
  119.     
  120.     node->sub.t = (Triangle *)calloc(1,sizeof(Triangle));
  121.     
  122.     node->sub.t->V[0] = &Q[ord(r1, r2, r3)];
  123.     node->sub.t->V[1] = &Q[ord(s1, s2, s3)];
  124.     node->sub.t->V[2] = &Q[ord(t1, t2, t3)];
  125.     
  126.     minmax(node->sub.t->V[0]->x,node->sub.t->V[1]->x,node->sub.t->V[2]->x,
  127.        &node->box[X][MIN],&node->box[X][MAX]);
  128.     minmax(node->sub.t->V[0]->y,node->sub.t->V[1]->y,node->sub.t->V[2]->y,
  129.        &node->box[Y][MIN],&node->box[Y][MAX]);
  130.     minmax(node->sub.t->V[0]->z,node->sub.t->V[1]->z,node->sub.t->V[2]->z,
  131.        &node->box[Z][MIN],&node->box[Z][MAX]);
  132.     
  133.     node->sub.t->pl = PlaneEquation(node->sub.t->V[0],
  134.                     node->sub.t->V[1],
  135.                     node->sub.t->V[2]);
  136.     
  137.   }
  138. }
  139.  
  140. /****************************************************************************/
  141.  
  142.  
  143. MeshSurf(m, d, level)
  144.      Mesh *m;
  145.      int d, level;
  146. {
  147.   Mesh *mp;
  148.   Edge *ep;
  149.   int i,j,k,l,I;
  150.   double r, s, tmp;
  151.   
  152.   if (m) {
  153.     if ((m->type == MESH) && (m->sub.m->type == FACE)) {
  154.       
  155.       mp = m->sub.m;
  156.       
  157.       if (mp->sub.e) {
  158.     
  159.     
  160.     if (level) {
  161.       
  162.       degree = d;
  163.       top = (degree+1)*(degree+2)*(degree+3)/6 - 1;
  164.       
  165.       for (i = degree; i > 0; i--)
  166.         for (j = degree-i; j >= 0; j--) {
  167.           k = degree-i-j;
  168.           
  169.           B[LOC(degree,i,j,k)] = *mp->sub.e->p;
  170.           B[LOC(degree,i-1,j+1,k)] = *mp->sub.e->next->p;
  171.           B[LOC(degree,i-1,j,k+1)] = *mp->sub.e->next->next->p;
  172.           
  173.           mp = mp->next;
  174.         }
  175.       
  176.       EDGE_DIVISIONS = 1 << level;
  177.       NUMBER_SAMPLES = (EDGE_DIVISIONS+1)*(EDGE_DIVISIONS+2)/2;
  178.       
  179.       Q = (Vertex *)calloc(NUMBER_SAMPLES,sizeof(Vertex));
  180.       
  181.       I = 0;
  182.       for (i = 0; i <= EDGE_DIVISIONS; i++) {
  183.         r = (double)i/EDGE_DIVISIONS;
  184.         for (j = 0; j <= (EDGE_DIVISIONS-i); j++) {
  185.           s = (double)j/EDGE_DIVISIONS;
  186.           
  187.           DECasteljau(r, s, 1-r-s);    
  188.           
  189.           for (l = top; l > top-4; l--)
  190.         for (k = 0; k < 4; k++)
  191.           B[l].P[k] /= B[l].P[W];
  192.           
  193.           Q[I].x = B[top].P[X];
  194.           Q[I].y = B[top].P[Y];
  195.           Q[I].z = B[top].P[Z];
  196.           
  197.           for (k = 0; k < 4; k++) {
  198.         B[top-1].P[k] = B[top-1].P[k] - B[top-3].P[k];
  199.         B[top-2].P[k] = B[top-2].P[k] - B[top-3].P[k];
  200.           }
  201.           
  202.           Q[I].a = B[top-1].P[Y]*B[top-2].P[Z] - B[top-2].P[Y]*B[top-1].P[Z];
  203.           Q[I].b = B[top-1].P[Z]*B[top-2].P[X] - B[top-2].P[Z]*B[top-1].P[X];
  204.           Q[I].c = B[top-1].P[X]*B[top-2].P[Y] - B[top-2].P[X]*B[top-1].P[Y];
  205.           
  206.           tmp = sqrt(Q[I].a*Q[I].a 
  207.              + Q[I].b*Q[I].b 
  208.              + Q[I].c*Q[I].c);
  209.           
  210.           Q[I].a = Q[I].a/tmp;
  211.           Q[I].b = Q[I].b/tmp;
  212.           Q[I].c = Q[I].c/tmp;
  213.           
  214.           I++;
  215.           
  216.         }
  217.       }
  218.       
  219.       Go(m, EDGE_DIVISIONS, 0, 0, 0, EDGE_DIVISIONS, 0, 0, 0, EDGE_DIVISIONS, level);
  220.       
  221.     }
  222.       }
  223.     } else if (m->type == MESH) {
  224.       
  225.       for (mp = m->sub.m; mp != (Mesh *) 0; mp = mp->next)
  226.     MeshSurf(mp, d, level);
  227.       
  228.       for (i = 0; i < 3; i++) {
  229.     m->box[i][MIN] = m->sub.m->box[i][MIN];
  230.     m->box[i][MAX] = m->sub.m->box[i][MAX];
  231.       }
  232.       
  233.       for (mp = m->sub.m->next; mp != (Mesh *) 0; mp = mp->next)
  234.     for (i = 0; i < 3; i++) {
  235.       if (mp->box[i][MIN] < m->box[i][MIN])
  236.         m->box[i][MIN] = mp->box[i][MIN];
  237.       if (mp->box[i][MAX] > m->box[i][MAX])
  238.         m->box[i][MAX] = mp->box[i][MAX];
  239.     }
  240.       
  241.     }
  242.   }
  243. }
  244.  
  245.  
  246.  
  247.